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Abstract 

The need for accurate material models to simulate the deformation, damage and failure of polymer 
matrix composites is becoming critical as these materials are gaining increased usage in the aerospace and 
automotive industries. While there are several composite material models currently available within LS- 
DYNA (Livermore Software Technology Corporation), there are several features that have been 
identified that could improve the predictive capability of a composite model. To address these needs, a 
combined plasticity and damage model suitable for use with both solid and shell elements is being 
developed and is being implemented into LS-DYNA as MAT_213. A key feature of the improved 
material model is the use of tabulated stress-strain data in a variety of coordinate directions to fully define 
the stress-strain response of the material. To date, the model development efforts have focused on 
creating the plasticity portion of the model. The Tsai-Wu composite failure model has been generalized 
and extended to a strain-hardening based orthotropic yield function with a nonassociative flow rule. The 
coefficients of the yield function, and the stresses to be used in both the yield function and the flow rule, 
are computed based on the input stress-strain curves using the effective plastic strain as the tracking 
variable. The coefficients in the flow rule are computed based on the obtained stress-strain data. The 
developed material model is suitable for implementation within LS-DYNA for use in analyzing the 
nonlinear response of polymer composites. 


Introduction 

As composite materials are gaining increasing use in aircraft components where impact resistance 
under high energy impact conditions is important (such as the turbine engine fan case), the need for 
accurate material models to simulate the deformation, damage and failure response of polymer matrix 
composites under impact conditions is becoming more critical. Within LS-DYNA (Ref 1), several 
material models are available for application to the analysis of composites. For example, the Chang- 
Chang failure model (Ref. 2) is utilized in MAT 22 and MAT 54. In these models, combinations of 
ratios of stresses to failure strengths are utilized to predict fiber or matrix based failure. The response is 
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assumed be linear elastic, with limited capability to simulate the nonlinear shear response. In MAT 22 
the failure is assumed to be brittle, while in MAT 54 the composite elastic constants are selectively 
reduced based on the failure mode, and a gradual unloading is permitted until ultimate element failure is 
reached. In MAT 58, a continuum damage model developed by Matzerrmiller et al (Ref 3) is employed, 
where the initiation and accumulation of damage is assumed to be the primary driver of nonlinearity in 
the composite response. The failure stresses and strains of the material in each of the coordinate directions 
are specified by the user, and the material stress-strain curves are approximated based on this data. The 
original version of the material model assumes a material response that is independent of strain rate. 
However, in the enhanced model MAT158, a modified version of the model used in MAT 58 is 
employed which incorporates strain rate dependence into the material response (Ref 1). A viscoelastic 
Prony series based on shear moduli is used to modify the computed stresses. The strain rate dependence in 
all of the coordinate directions is assumed to be identical. MAT161 is an additional continuum damage 
mechanics model available within LS-DYNA (Ref. 4). In this model, a stress to strength ratio approach in 
a variety of coordinate directions is used to specify the beginning of either fiber or matrix based failure. 
Appropriate failure mode based damage functions are then used to compute the reduction in elastic 
moduli in each of the coordinate directions. In MAT 219, based on the CODAM model developed at the 
University of British Columbia (Ref 5), a strain versus failure strain ratio based approach is used to 
predict the initiation of damage in a sublaminate of the composite. Separate damage accumulation and 
modulus reduction functions based on failure mode and coordinate direction are implemented in a 
continuum damage mechanics formulation. In MAT 22 1 , damage accumulation functions based on 
current strains, damage initiation strains and failure strains are used to reduce the elastic moduli of the 
composite in each of the coordinate directions in a damage mechanics approach (Ref 1). In MAT 26I, 
which is based on a damage and failure model developed by Pinho et al (Refs. 6 and 7), separate models 
for fiber tension failure, fiber kinking failure, matrix tensile failure and matrix compressive failure are 
developed, where different functional forms are used for each of the failure criteria (for example, the 
matrix compression failure criterion is based on an extension of the Mohr-Coulomb failure criterion). The 
various failure criteria are combined together using fracture mechanics concepts to develop a constitutive 
model. In MAT 262, based on a model developed by Camonho, et al (Refs. 8 and 9), an energy approach 
is utilized to generate damage functions in various coordinate directions within the context of a 
continuum damage mechanics formulation. 

While all of the material models discussed above have been utilized with some level of success in 
modeling the impact response of polymer composites, there are some areas where the predictive 
capability can be improved. In general, the existing models often require significant a priori knowledge of 
the damage and failure such that their use as predictive tools can be limited. While these models generally 
assume that the nonlinear response of the composite is due to either deformation mechanisms (plasticity) 
or damage mechanisms, an improved model should have the capability to simulate the actual material 
behavior in which the material nonlinearity is due to a combination of both deformation and damage 
mechanisms. The input to the current material models generally consists of point-wise properties that lead 
to curve fit approximations to the material stress-strain curves. This type of approach leads to either 
models with only a few parameters, which provide a cmde approximation at best to the actual stress-strain 
curve, or to models with many parameters which require a large number of complex tests to characterize. 
An improved approach would be to use tabulated data from a well-defined set of experiments to 
accurately define the complete stress-strain response of the material. Furthermore, many of the existing 
models are only suitable for use with two-dimensional shell elements, which cannot capture the through- 
thickness response which may be significant in impact applications. Ideally, a fully three-dimensional 
formulation suitable for use with solid elements would be available, along with the shell element 
formulation. In addition, incorporating the effects of strain rate is only possible in a few of the existing 
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models, and in those models, the effects of strain rate on the material response are assumed to be the same 
in all of the material coordinate directions. 

To begin to address these needs, a multi-institution consortium has been formed to develop and 
implement a new composite material model within LS-DYNA, which will be implemented as MAT 213. 
Currently, the primary focus of the effort has been on the development and testing of the orthotropic 
macroscopic three-dimensional plasticity model. As will be discussed in more detail in the following 
sections of this paper, the commonly used Tsai-Wu composite failure criterion (Ref 10) has been 
generalized and extended to a strain-hardening model with a quadratic yield function and a nonassociative 
flow rule. A similar approach has been used by Sun and Chen (Ref 1 1), in which a general quadratic plastic 
potential function was generalized into a plasticity model suitable for use with composites. The coefficients 
of the yield function for the new composite model are determined based on tabulated stress-strain curves in 
the various normal and shear directions, along with selected off-axis curves. The non-associative flow rule is 
used to compute the components of the plastic strain along with the effective plastic strain. The stresses that 
are used in the flow mle are computed based on the input tabulated stress-strain curves. Systematic 
procedures have been developed to compute the various coefficients in the yield function and the flow rule 
based on the tabulated input data. An important point to note is that the developed material law is not limited 
to the analysis of unidirectional plies alone. Alternatively, the material model is meant to be a fully 
generalized model suitable for use with any composite architecture (unidirectional, laminated or textile). In 
the following sections of this paper, a detailed derivation of the yield function and the flow law are 
presented. The procedures to be used to characterize the material constants in the yield function and the flow 
law are also discussed in detail. Details of the numerical implementation of the material model within LS- 
DYNA and descriptions of verification and validation studies that were conducted using the developed 
material model are available in a companion paper (Ref 12). 


Material Law Derivation 

A general quadratic three-dimensional orthotropic yield function based on the Tsai-Wu failure model 
is specified as follows, where 1, 2, and 3 refer to the principal material directions. 
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In the yield function, o^ represents the stresses and Fjj and are coefficients that vary based on the 
current values of the yield stresses in the various coordinate directions. By allowing the coefficients to 
vary, the yield surface evolution and hardening in each of the material directions can be precisely defined. 
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The values of the normal and shear coefficients can be determined by simplifying the yield function for 
the case of unidirectional tensile and compressive loading in each of the coordinate directions along with 
shear tests in each of the shear directions, with results as shown below. 
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In the above equation, the stresses are the current value of the yield stresses in the normal and shear 
directions (determined using procedures to be discussed below), where the superscript T indicates the 
tensile yield stress, and the superscript C indicates the absolute value of the compressive yield stress. To 
determine the values of the off-axis coefficients (which are required to capture the stress interaction 
effects), the results from 45° off-axis tests in the various coordinate directions can be used. For example, 
the stresses in the local material axis system resulting from a uniaxial tensile test of a [45°] composite (for 
a unidirectional composite), or any uniaxial tensile test conducted at 45° in the 1-2 plane from the 
longitudinal (1) material axis for a multiply laminated or textile composite, can be computed to be the 
following by the stress transformation equations (Ref 10). 


an =0.5ay45 

G 22 = 0.5ar45 (3) 

ai 2 = -0.5ay45 


In the above equation, 0^45 is the structural level current uniaxial yield stress resulting from a tensile test 
on a [45°] ply, and the remaining stresses are the material axis system based stresses. By substituting this 
equation into Equation (1) and solving for the off-axis constant Fj2, the following expression can be 
obtained 


Fn^ . ^ -Ufh+Fii +F44) (4) 

(ay45j ^ 

where all of the terms are as defined earlier. By using similar procedures involving 45° rotations in the 2- 
3 plane and the 1 -3 plane, similar expressions for the constants and can be determined. 

A nonassociative flow rule is used to compute the evolution of the components of plastic strain. The 
plastic potential for the flow rule is shown below 

h = 

I (5) 

+H22<J22 +'^33a|3 +2Hi 2<JU<J22 +2H23<J22<J33 +2H3\0330n +H 440 I 2 + FI 550^2 +2f66a|j 

where Oy are the current values of the stresses and //y are independent coefficients, which are assumed to 
remain constant. Procedures for determining the values of the coefficients will be discussed later in this 
paper. The plastic potential function in Equation (5) is used as follows in a flow law, where the usual 
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normality hypothesis from classical plasticity (Ref. 13) is assumed to apply and the variable 1 is a scalar 
plastic multiplier. 
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In the above equation, kfj are the components of plastic strain rate. For the shear components, tensorial, 

not engineering, strains are used. By examining ratios of the components of plastic strain in response to a 
given uniaxial stress component, the normal and off-axis coefficients 77^ in the plastic potential function 

can be related to “plastic Poisson’s ratios”, v^ , as follows. This concept, as will be shown later in this 

paper, is useful in developing procedures to characterize the values of the coefficients. 


011 



022 ^ 0 



033 ^ 0 


^22 

Hu 

yP = 

pP 

Hu 

yP 

II 

•CO 1 
1 

II 

H22 

^11 

Hn 

21 

pP 

t-22 

H 22 

'22 

pp 

^33 

Hu 

pP 

^33 

Hu 

V?. - 

pP 

^33 _ 

7/23 


pp 

^11 _ 

Hu 

pP 

^11 

Hn 

'^23 

pP 

^22 

H 22 

'^31 

pp 

^33 

H33 


(7) 


In the above equation the expressions in the first column are only valid for the case where the only 
nonzero stress is in the 1 direction, the expressions in the second column are only valid for the case where 
the only nonzero stress is in the 2 direction, and the expressions in the third column are only valid for the 
case where the only nonzero stress in in the 3 direction. An important point to note is that while the yield 
function can accommodate differences between the tensile and compressive responses of the composite, 
the flow law, due to its lack of linear terms, cannot make such a distinction. However, introduction of a 
linear term would make the plastic Poisson’s ratio, and thus the coefficients in the flow law, highly 
dependent on the stress which would most likely lead to erratic behavior of the model. Another 
consideration is that the coefficients of the flow law could be assumed to vary based on the current stress 
and strain state, but that would require developing evolution conditions on the flow law, which would be 
very difficult to determine based on experimental data. 

Given the flow law, the principal of the equivalence of plastic work (Ref 1 3) can be used to 
determine expressions for the effective stress and effective plastic strain. By taking the vector product of 
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the stress tensor and the plastic strain rate tensor, one ends up with the plastic potential function 
multiplied by the plastic multiplier k . Since, by the principal of the equivalence of plastic work, the 
vector product of the stress tensor and the plastic strain rate tensor must equal the product of the effective 
stress and the effective plastic strain rate, one can conclude that the plastic potential function h can be 
defined as the effective stress and the plastic multiplier k can be defined as the effective plastic strain 
rate. This concept is expressed mathematically in the equation shown below 

• Bh 

WP -a:i.P -a:k — -hk-<3e&e (8) 

5o 

where Og is the effective stress and is the effective plastic strain rate. 

To compute the current value of the yield stresses needed for the yield function, the common practice 
in plasticity constitutive equations is to use analytical functions to define the evolution of the stresses as a 
function of the components of plastic strain (or the effective plastic strain). Alternatively, in the 
developed model tabulated stress-strain curves are used to track the yield stress evolution. The user is 
required to input twelve stress versus plastic strain curves. Specifically, the required curves include 
uniaxial tension curves in each of the normal directions (1,2,3), uniaxial compression curves in each of 
the normal directions (1,2,3), shear stress-strain curves in each of the shear directions (1-2, 2-3 and 3-1), 
and 45° off-axis tension curves in each of the 1-2, 2-3 and 3-1 planes. The 45° curves are required in 
order to properly capture the stress interaction effects. By utilizing tabulated stress-strain curves to track 
the evolution of the deformation response, the experimental stress-strain response of the material can be 
captured exactly without any curve fit approximations. 

The required stress-strain data can be obtained either from actual experimental test results or by 
appropriate numerical experiments utilizing stand-alone codes. Currently, only static test data is 
considered. Future efforts will involve adding strain rate and temperature dependent effects to the 
computations. For certain composite architectures, the number of stress-strain curves that need to be input 
can be simplified. For example, for unidirectional composites, the normal response in the 3 direction can 
be assumed to be identical to the normal response in the 2 direction due to the transverse isotropy of the 
material. Likewise, the shear response in the 1-3 direction can be assumed to be identical to the shear 
response in the 1-2 direction. Similarly, the 45° off-axis curve in the 1-3 plane can be assumed to be 
identical to the 45° off-axis curve in the 1-2 plane. Due to the transverse isotropy of a unidirectional 
composite, one can also demonstrate that for the case of a unidirectional composite the 45° off-axis curve 
in the 2-3 plane is approximately equal to the normal stress-strain curve in the 2 (or 3) direction. For other 
composite architectures with significant symmetry (such as two-dimensional plain or satin weave 
composites), a certain degree of simplification in the stress-strain data that is required to be input is also 
possible. 

To track the evolution of the deformation response along each of the stress-strain curves, the effective 
plastic strain is chosen to be the tracking parameter. Using the numerical procedure described briefly later 
in this paper and in more detail in a companion paper (Ref 12), the effective plastic strain is computed for 
each time/load step. The stresses for each of the tabulated input curves corresponding to the current value 
of the effective plastic strain are then used to compute the yield function coefficients. To enable this 
tracking, each of the input stress-strain curves needs to be converted into a plot of stress versus effective 
plastic strain. To enable this conversion, the principal of the equivalence of plastic work is used again to 
relate the unidirectional plastic strain increment to the effective plastic strain increment. For example, for 
the case of unidirectional loading in the 1 direction the effective plastic strain increment (and integrated 
current total effective plastic strain) can be computed as follows 
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where Ou is the unidirectional stress in the loading direction, de[^ is the plastic strain increment in the 

loading direction and dSg is the increment in effective plastic strain. Note that in a general case the 
increment in effective plastic strain can be computed as the incremental area under the stress-plastic strain 
curve divided by the current value of the effective stress. A graphical representation of this process is 
shown below in Figure 1 . 

To compute the evolution of the effective plastic strain, a numerical algorithm based on the radial 
return method is employed. Given the revised value of the effective plastic strain, the yield stress values 
and the overall stress state of the material can be updated. Details of the numerical implementation are 
provided in a companion paper (Ref 12), but some of the key theoretical fundamentals of the algorithm 
are given below. The standard elastic constitutive equation is used to compute the revised stresses, where 
the flow law is applied for the computation of the plastic strains. 


o = C:(e-£p) = C: 
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( 10 ) 


In the above equation, C is the standard elastic stiffness matrix, e is the total strain rate tensor, and the 
remaining terms are as defined previously. To compute the effective plastic strain rate, the consistency 
condition is applied in combination with the elastic constitutive equation 



User provided load curves are Internally stored is true stress 

true stress versus true plastic strain versus effective plastic strain 

Figure 1. — Conversion of stress versus plastic strain curves to stress versus effective 
plastic strain curves. 
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where /is the yield function defined in Equation ( 1 ) and q is the vector of yield stresses based on the 
input stress-strain curves. 

q = (afj 022 *^33 ‘^n *^22 *^^3 *^12 cr 23 031 0745 0745-23 0745-31) ( 12 ) 


where 0745-23 is the yield stress resulting from a uniaxial tensile test conducted at 45 ° in the 2-3 plane 
and 0745-31 is the yield stress resulting from a uniaxial tensile test conducted at 45 ° in the 3-1 plane. To 
compute the derivative of the components of the q vector with respect to the effective plastic strain A,, a 
chain rule type of computation is carried out where the derivatives of the individual components of the 
vector (the respective yield stresses) with respect to the effective plastic strain are determined as follows. 


dk dzfj dk 

The derivative of the components of the yield stress vector with respect to plastic strain is the 
instantaneous slope of the stress versus plastic strain curve from each of the 12 input stress-strain curves. 
The derivative of the plastic strain with respect to the plastic multiplier k can be computed from the flow 
law for each of the specialized unidirectional (or off-axis) cases represented by the input stress-strain 
curves. Note that for the off-axis cases the flow law must first be written in the material axis system and 
then the plastic strains must be transformed into the structural axis system in order to compute the proper 
derivative of the plastic strain with respect to the plastic multiplier. A secant iteration approach is used in 
combination with the radial return method to compute the specific value of the effective plastic strain rate 
(Ref 12 ). 


Characterization of Flow Law Coefficients 

The coefficients in the flow rule need to be characterized based on data obtained from experimental 
stress-strain curves for the composite. Specifically, data from the 12 experiments leading to the input 
stress-strain curves must be utilized. The data can be obtained either from actual experiments or 
numerical experiments. For example, if one knows the mechanical properties and stress-strain curves for 
the fiber and matrix constituents, numerical analyses can be conducted using either high fidelity finite 
element analysis (where the fiber and matrix are modeled as distinct constituents) or analytical tools such 
as the NASA Glenn developed micromechanics code MAC/GMC (Ref 14 ) to obtain numerical 
approximations to the actual composite experimental data. 

First, the procedures used to obtain the required coefficients for the simplified case of a unidirectional 
carbon fiber based polymer matrix composite will be discussed. Afterward, a generalized procedure to 
obtain the required constants for a general composite will be presented. For a unidirectional carbon fiber 
composite, as discussed by Sun and Chen (Ref 1 1 ) a reasonable approximation to make is that due to the 
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linear behavior of the carbon fiber the plastic strain in the fiber direction (the 1 direction) is equal to zero 
for all values of stress. By examining the second expression in Equation (6), one can conclude that for a 
general stress state the only way the plastic strain in the 1 direction can be zero is for the coefficients 
Hi2 and to all be identically zero. Next, since the transverse tension response of the composite can 
display a certain degree of nonlinearity, a reasonable approximation to make is that for the case of a 
unidirectional load in the 2 direction the effective stress is equal to the applied transverse tension stress 
O22. By simplifying the plastic potential function h (shown in Equation ( 5 )) for the case of a unidirectional 
transverse tension stress, and by setting this term (along with the transverse tensile stress) equal to the 
effective stress, one can show that the value of the coefficient H22 must equal one. Since a unidirectional 
composite is transversely isotropic, one can set the value of to also equal one. Given these values, 

using the definition of the plastic Poisson’s ratio v^3 shown in Equation ( 7 ), one can find that the 

coefficient H22 equals the negative of the plastic Poisson’s ratio v^3 . Since the flow surface coefficients 

are assumed to be constant, an average value of the plastic Poisson’s ratio must be determined from the 
experimental and/or numerical test data from a unidirectional transverse tension test. To determine the in- 
plane shear coefficient 7/44, utilizing a procedure similar to that developed by Sun and Chen (Ref 1 1 ) for 
their model, the plastic potential function h given in Equation ( 5 ) and the plastic strain definition given in 
Equation (6) can be simplified for the case of in-plane shear loading in the 1-2 direction, yielding the 
following expressions. 


h-<3e- 

, pP ( 14 ) 

sf2 = O.Ssf ^ 

Starting from the experimental stress versus plastic strain plot for the case of in-plane shear loading, 
using the expressions in Equation ( 14 ) the value of 7/44 can be optimized to provide the closest possible 
match with the overall effective stress versus effective plastic strain curve, selected for the case of a 
unidirectional carbon fiber composite to be the transverse tension stress versus plastic strain curve. Due to 
the transverse isotropy of the unidirectional composite, the value of the coefficient Tfgg can be set equal to 
the value of the coefficient 7/44. To determine the value of the coefficient, a similar procedure can be 

carried out for the case of shear loading in the 2-3 direction. Alternatively, due to the transverse isotropy 
of the composite, by computing the effective stress for the case of off-axis loading in the 2-3 plane, an 
isotropic type of relation can be determined for the coefficient 

7755= 2(1 + v^3) ( 15 ) 

where v^3 is the plastic Poisson’s ratio in the 2-3 direction, obtained from the material data. 

Similar procedures are used for the case of a general composite with no identified symmetries. An 
example of a general composite would be a triaxially braided composite which is modeled as a smeared, 
homogenous material. Another example would be a generalized laminated composite modeled as a 
smeared, homogeneous material. For the case of a general composite, a unidirectional tension test in the 
longitudinal ( 1 ) direction is chosen as the baseline case. By utilizing Equation ( 5 ) and setting the plastic 
potential function and the unidirectional longitudinal stress equal to the effective stress, the value of the 
77 ii coefficient can be found to equal one. Based on the relations given in Equation ( 7 ), the value of the 

coefficient H12 can be found to equal the negative of the plastic Poisson’s ratio v^2 > the value of the 
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coefficient //j3 can be found to be equal to the negative of the plastic Poisson’s ratio . As before, the 

value of the plastic Poisson’s ratio used in the computations is the average value of the parameter 
determined from the experimental and/or numerically obtained data. In a similar fashion, the values of the 
coefficients H 22 , /^23 can be computed. 


H22 = 


'12 


*'21 


//23=-vf3^ 


'21 


//33 = — 
P 

V3I 


(16) 


To compute the values of the coefficients I/ 44 , and H(^^, the same procedure that was utilized to 

compute H 44 for the case of a unidirectional carbon fiber composite is used, only all three shear loading 
cases (012 , c^ 23> 031) must be used. 


Conclusions 

A generalized plasticity model suitable for use in polymer composite impact simulations has been 
developed. The plasticity model will be part of a generalized plasticity and damage model which will be 
implemented into LS-DYNA as MAT 213. The yield function is based on the Tsai-Wu composite failure 
model, and a suitable nonassociative flow rule was defined. Methods of utilizing tabulated stress-strain 
data to track the evolution of the yield stresses as a function of the effective plastic strain have been 
developed. The elastic constitutive equation can be used in combination with the consistency condition to 
compute the effective plastic strain. In a companion paper (Ref 12), further details of the numerical 
algorithm are provided along with examples of verification and validation studies which have been 
performed to examine the accuracy and efficiency of the analytical model. 

Future efforts will involve the development of an equivalent damage model, to be employed in 
combination with the plasticity model, which will also be based on utilizing tabulated stress-strain data to 
track the evolution of damage in the composite. Element failure will also be incorporated within the 
model. In addition, strain rate and temperature effects will be added to the analysis capability. Overall, 
when completed the composite model as implemented in MAT 213 will provide significant 
improvements to the state of the art in the modeling of the impact response of polymer matrix composites. 
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